Nano-engineered non-uniform strain in graphene 
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Recent experiments showed that non-uniform strain can be produced by depositing graphene over 
pillars. We employed atomistic calculations to study the non-uniform strain and the induced pseudo- 
magnetic field up to 5000 Tesla in graphene on top of nano-pillars. By decreasing the distance 
between the nano-pillars a complex distribution for the pseudo-magnetic field can be generated. 
Furthermore, we performed tight-binding calculations of the local density of states (LDOS) by 
using the relaxed graphene configuration obtained from the atomistic calculations. We find that the 
quasiparticle LDOS are strongly modified near the pillars, both at low energies showing sub-lattice 
polarization, and at high energies showing shifts of the van Hove singularity. Our study shows 
that changing the specific pattern of the nano-pillars allows us to create a desired shape of the 
pseudo-magnetic field profile while the LDOS maps provide an input for experimental verifications 
by scanning tunneling microscopy. 



Graphene is a newly discovered atomic thin two- 
dimensional honeycomb lattice consisting of carbon 
atomsi. It is a zero gap semimetal with a conical band 
structure where the conduction and valence bands touch 
each other at the Dirac point^. Nano-engineered non- 
uniform strain distribution in graphene is a promising 
road to generate a band gap and a pseudo-magnetic 
field. Scanning tunneling microscopy (STM) measure- 
ments have shown strain-induced Landau levels^ which 
correspond to a large pseudo-magnetic field. Shear strain 
is essential and neither uniaxial nor isotropic strain pro- 
duces a strong uniform pscudo- magnetic field 3 . 

Graphene's highly responses to external forces result- 
ing in mechanical deformations. Over the last few years 
there have been many efforts to control graphene's elec- 
tronic properties by strain^—. Elastic deformations cre- 
ate a pseudo-magnetic field which acts on graphene's 
massless charge carriers^—. The resulting variation 
of the hopping energies can be viewed as an induced 
pseudo-magnetic field which enters in the Dirac equa- 
tion. Engineering of the right topology of the induced 
pseudo-magnetic field provides symmetrical magnetic 
confinement which confines electrons in specific regions 
in spaced. 

Recently, it was predicted that non-uniform strain may 
lead to a considerable energy gap and a large gauge field 
that effectively acts as a uniform magnetic field — . Re- 
cently, Tomori et al used pillars made of a dielectric ma- 
terial (electron beam resist) which were placed on top 
of a substrate which is then overlayed with graphene 
to generate non-uniform strain on a micro-scale^ 3 -. The 
graphene sections which are located between the pillars 
are attached to the substrate and the size and separa- 
tion of the pillars control the strength and distribution 
of the strain. The length scale in the experiment was 
micronmctcres and Si02 was used as the substrate^ 3 -. 

Here we study non-uniform strain at the atomistic 
scale where the continuum approach is no longer appli- 
cable. We also study the local density of state (LDOS) 
maps using the relaxed graphene configuration as input 



for tight-binding calculations. We find very strong non- 
uniform pseudo-magnetic fields that can be created by 
depositing graphene on a substrate decorated with nano- 
scale pillars and find that the quasiparticle LDOS are 
strongly modified near the pillars. The optimum config- 
uration of graphene over such nano-pillars depends on 
the imposed boundary conditions. The induced pseudo- 
magnetic fields are larger than 1000 Tesla and are spa- 
tially distributed around the nano-pillars. Decreasing the 
distance between the nano-pillars alters the six-fold sym- 
metry of the pseudo-magnetic field distribution and re- 
sults in a new configuration of magnetic confinement for 
the charge carriers on the graphene around the nano- 
pillars. Our study shows the LDOS maps around the 
pillars, which can be experimentally verified by STM. 

Atomistic model. Classical atomistic molecular dy- 
namics simulation (MD) is employed to find the optimum 
configuration of large flakes of graphene (GE) over the 
nano-pillars. The second generation of Brenner's bond- 
order potential is employed and is able to describe cova- 
lent sp 3 bond breaking and the formation of associated 
changes in atomic hybridization within a classical poten- 
tial^. The van der Waals (vdW) interaction between GE 
and nano-pillars/substrate is modeled by employing the 
Lennard- Jones (LJ) potential^—. 

In order to model the substrate, a (100) surface with 
lattice parameter equal to t=2> A is assumed with L J pa- 
rameters as and 65. The density of the sites in the sub- 
strate is S5 = £~ 2 and the number of atoms is 13700. 
Nano-pillars are double- wall armchair carbon nanotubes 
(DWCNT) taken with (3,3) and (6,6) indexes including 
144 atoms (see left insets in Fig[T]). The number of atoms 
in the graphene sheet is 44800 which is equivalent to a 
sheet of size 34.8 x 34.43 nm 2 . We assume that both the 
substrate and nano-pillar atoms are rigid during the sim- 
ulation. 

To model the interaction between two different types 
of atoms such as carbon atom (C) and substrate atom 
(S), we adjust the LJ parameters using the equations 
e T = \/£c e s an d or = (o'c + cr s)/2. The parameters for 



carbon arc gq — 

3.369 A and e c = 2.63 meV. For the 
substrate atoms we took 175 =3.5A and £5 =10.0 meV 
which are typical for insulators, e.g. SiC>2^. 

The atomic stress experienced by each i th atom can be 
expressed a a 19 i 20 

where the inner summation is over all the carbon atoms 
which are neighbors of the i th atom which occupies a 
volume O = Attciq/S. The quantities m and v l denote 
the mass and velocity of i th atom and the scaler r"j is 
the v component of the distance between atoms 'i' and 
T- Ftj i s the force on i th atom due to atom j th in the /1 
direction. We used this expression to calculate the stress 
on each atom. In order to be able to visualize the stress 
distribution on the GE atoms, we colored the atoms using 
a dimensionless invariant quantity^ J2 = \[{Vxx~ Vyy) 2 ~ s r 

(Vyy ~ Vzz) 2 + {Vzz - Vxx) 2 + 6{vly + viz + Viz)] > i- e - dark 
green (white) is related to a minimum (maximum) value 
of J 2 . 

Strain induced pseudo-magnetic field. Coupling 
the Dirac equation, which governs the low energy elec- 
tronics of graphene, to the curved surface is a common 
way to study the effects of graphene's curved geometries 
on its corresponding electronic propertie a 9 ' 10 . The metric 
of the curved surface describes the curvature of the sur- 
face. The origin of the deformations are external stresses 
that deform graphene so that the nearest neighbor dis- 
tances become non-equal. The latter results in modi- 
fied hopping parameters, which are now a function of the 
atomic positions £(r)ii. Assuming small atomic displace- 
ments (u = rj— Ti < ao) and rewriting the Dirac Hamilto- 
nian in the effective mass approximation with non-equal 
hopping parameters yields the induced gauge fields: 

_ 2/3H 

A- — \U X x ^yyi 2li;ry), (2) 

oaoe 

where ~ 3 is a constant and u a p is the strain tensor in- 
cluding out of plane displacements 9 -. The corresponding 
pseudo-magnetic field perpendicular to the x — y plane is 
obtained as 

B = dyA x -d X Ay. (3) 

This is the pseudo-magnetic field that an electron expe- 
riences in the K-valley. We will find B by making the 
necessary differentiations numerically in the case of sup- 
ported boundary conditions. Here we found that the ma- 
jor contribution is due to the out-of-planc terms, which 
appear mainly around the deformed parts. The other in- 
plane terms contribute less to the pseudo-magnetic field 
around the deformed parts particularly when the system 
is large as compared to the deformed regions. 

Tight binding model .The electronic properties are 
described by a tight-binding Hamiltonian for the 7r car- 
bon orbitals. The minimal Hamiltonian, which describes 
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FIG. 1: (Color online) The optimal configuration of graphene 
on top of double wall armchair carbon nanotubes (nano- 
pillars) and a square lattice substrate (right-insets shows a 
top view). The left-insets show the nano-pillars and the sub- 
strate. In (a) the distance between pillars is 10 nm while in 
(c-d) it is 5 nm. All pillars are at the same height except (c) 
where the central pillar is 2nm higher. Colors indicate the 
scaled stress distribution, i.e. white represents highest stress 
and dark-green lowest stress. 



the low-energy band structure is: 

n = E -KnA&Cj,, + h.C (4) 

where ct (cj a ) creates (destroys) an electron at site % 
(J). The sum runs over nearest neighbors pertaining to 
opposite sub-lattices (i,j) and the electron spin, a. In 
the following we will ignore the spin degrees of freedom 
since no spin-flipping term is present in the Hamiltonian. 

The strain is included in the modified hopping ampli- 
tudes between 7r orbitals, t w (rij), according to the empiri- 
cal relation i^(ry) = 70 cxp 3 ' 37 '~ 1 \ where 70 = 2.7eV 
and ao = 1.42A is the equilibrium inter-carbon distance^. 
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FIG. 2: (Color Online), (a-d) Pseudo-magnetic fields in one Dirac cone for 4 different pillar configurations shown in Fig. 1, 
respectively, (e-h) Low energy (E=0.1332eV) LDOS map and (i-1) high energy (E=2.2306eV) LDOS map for the same 4 
configurations. Note that the bulk unstrained LDOS is subtracted from the LDOS maps. 



This also gives a good approximation for the next-nearest 
neighbor hopping amplitude. We also consider the effect 
of misalignment of the tt orbitals due to the finite curva- 
ture. This effect translates into the mixing of the tt and a 
orbitals. Depending on the local curvature, the modified 
hopping amplitude is: 

t(rij) = t n sm(8i) sm(8j)cos(4>) — t a cos(#i) cos(#j), (5) 

where 6i and 6j are the angles formed by the normals at 
each atom, n, and rij, with the inter-atomic distance 
and <p describes the angle formed by the normal rij and 
the plane defined by and r^^Sr— . The strain configu- 
ration and the curvature are extracted from the relaxed 
position of the graphene sheet obtained from our molec- 
ular dynamics simulation. For the systems considered 
here the effect of curvature is small when compared to 
the effect of strain on the hopping amplitudes. 

Since the strain is inhomogeneous it is not possible to 
use any symmetries in the calculations of the electronic 
properties. Therefore, the system size considered here 
(44800 carbon atoms) becomes prohibitively large for an 



exact digitalization of the Hamiltonian. Instead we will 
numerically obtain an approximation of the Green's func- 
tion by using a Chebyshev expansion within the Kernel 
Polynomial method^ - — . The Green's function is defined 
as: 

Gijiu) = (x\d(u,)\4) (6) 

where G(u + iff) = [uj + it] — 

First a scaling of the excitation energies is performed, 
e.g. H = {H — lb) /a, Q = (u> — b)/a where a = (E max — 
E m i n )/(2 - rj) and b = (E max + E mm )/2 where r) > 
is a small number. Following Refs. [2al27j . the Green's 
function's components can be expressed as an expansion 
written in terms of Chebyshev polynomials: 

_„. oo 

Gyp) = -_= £ a n (i, j)e ™ (*> (7) 

where the coefficients a n (i,j) = (ci\v n ) can be obtained 
by an iterative procedure involving repeated applications 
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of the Hamiltonian on iterative vectors \v n ): 

\v n +i) = 2n\v n ) - (8) 

where |uo) = |cj) and \v-%) = 0. Significant computa- 
tional speed-up is achieved when the computations are 
done on graphical processing units (GPU), i.e. video 
cards. The computations are performed on Nvidia 
GeForce GTX 580 cards. 

The physical properties that can be straightforwardly 
extracted from the Green's functions are the local density 
of states (LDOS): 

Ni(u) = --/m[Gf«(w)], (9) 

the factor of 2 appears due to the summation over the 
spin components. 

Results and discussion. At the start of our sim- 
ulation we put graphene on top of the nano-pillars at 
ho = 1.4 nm. The substrate is at zero height and the 
nano-pillars arc located between graphene and the sub- 
strate. We have investigated two particular patterns of 
nano-pillars: i) five DWCNTs which have in-plane co- 
ordinates (0,0) and (±d,±d) with d = 5,10nm, and ii) 
four DWCNTs at (0,0), (±dV3/2,±d/2) and (0,-d). The 
height of DWCNTs was set to be 1 nm (except in Fig.[TJc) 
which is 1.5 nm). In order to prevent crumpling at the 
boundaries we only allow the boundary atoms to vibrate 
in the z-direction. 

In Fig. [T] we show the optimal configuration of GE on 
top of four (a-c) and five (d) DWCNTs. Right insets in 
Figs. nja-d) show a top view. Notice that the stress dis- 
tribution is mainly concentrated around the nano-pillars 
as expected. For the configuration presented in Fig. Ufa) 
the pillars are lOnm apart. Due to the vdW interac- 
tion the graphene sheet will stick to the substrate except 
around the pillars where the shape is close to a Gaus- 
sian even though deviations from an isotropic descrip- 
tion exists. In the atomic limit, a slight anisotropy ap- 
pears, the graphene sheet bends mostly in the zig-zag 
direction making the shape of the deformation hexago- 
nal. In Figs. [TJb-d) the pillars are closer together (i.e. 
5nm). Due to its large bending rigidity, the graphene 
sheet will be suspended over the substrate in the regions 
between the pillars. Depending on the pillar configura- 
tion, various stain configurations are achieved. If all the 
pillars have the same height, Figs. [IJb) and 1(d), most 
of the stress is obtained at the pillar location and where 
graphene sticks to the substrate. If one of the pillars is 
higher, Figs. HJc), besides the maximal stress at the pil- 
lar location, high stress is also obtained throughout the 
suspended sheet. 

The corresponding pscudo magnetic field profiles gen- 
erated by the strain configurations are shown in Figs.^a- 
d). When the deformations are isolated, Fig. HJa), 
the gauge field and the pseudo magnetic field exhibit 
six fold symmetry&ii similar to a Gaussian deforma- 

2, 2 

tion, h(x,y) = Gcxp(— X 2 ^ ). The continuum the- 
ory predicts that the pseudo-magnetic gauge field is 



A = h ^£i (x 2 — y 2 > ~ x u) an d the pseudo-magnetic field 

is B = fefe^ (x 2 + y 2 ) sin(36'), where 6 is the azimuthal 
angle. Large pseudo-magnetic fields on the order of thou- 
sands of Teslas are obtained. When the graphene sheet is 
suspended, Figs. HJb-d), the six-fold symmetry survives 
near the pillars but more complex pseudo-magnetic field 
profiles can be obtained, from large fields throughout the 
suspended sheet, Fig. GJc), to fields localized only near 
the edges of the suspended sheet, Fig.[2Jd). As seen from 
Figs. [2jc,d) the closer the pillars, triangular and rectan- 
gular magnetic field profile is created within the position 
of the pillars. In Figs. &c) there is a high magnetic field 
region at the center and the electron can not pass through 
this region. 

In order to investigate the effect of the strain on the 
electronic properties, we input the relaxed positions of 
the atoms obtained from the atomistic simulation into 
the tight-binding model in order to find the LDOS maps 
around the pillars. These are shown in Figs. [UJe-h) for 
E = 0.1332eV and in Figs. C^i-l) for E = 2.2306eV. 
Two regimes can be observed, depending on the en- 
ergy. For low energies, the pseudo-magnetic field will 
induce sub-lattice polarized states localized either near 
the pillars, Fig.^e), or in the regions with large pseudo- 
magnetic fields, Fig. [2jg). In the five-pillar configura- 
tion, Fig. [UJh), these low energy states are mostly local- 
ized near the edges of the suspended region. Interference 
patterns which depend on the energy are observed^. A 
very different effect, which is not described by the low 
energy Dirac approximation^, is related to the shift of 
the van Hove singularity seen in unstrained graphene at 
E = 70 = 2.7eV. Because of strain, the hopping pa- 
rameters will be modified, therefore locally shifting the 
van Hove singularity. This is observed in Figs. HJi-l), 
where the enhancement of the LDOS is correlated with 
the stress and is enhanced where the stress is larger. Ad- 
ditional interference patterns appear between the pillars. 
The deviation from the Gaussian shape of the LDOS 
modification, i.e. the hexagonal shape, of the isolated 
deformations is also obvious from Fig. [2ji) . 

Conclusions. By combining molecular dynamics 
simulations with tight binding calculations we have 
shown how strain can be manipulated at the nano- 
scale. Isolated pillars show a six-fold symmetric pseudo- 
magnetic field and LDOS modification. By decreasing 
the distance between the pillars, the six fold symmetry 
of the pseudo-magnetic field is altered and a complex 
field profile appears within the suspended area. We found 
that by modifying the inter-pillar distances and the pillar 
heights one can design a particular desired magnetic field 
profile. Modifications of the hopping parameters due to 
changes in the C-C distances induced by strain, mod- 
ify the LDOS around the deformations of the graphene 
sheet. Verifications of the six-fold symmetry of the LDOS 
near pillars could be easily confirmed with STM experi- 
ments. 

Acknowledgment. This work was supported by the 
Flemish Science Foundation (FWO-V1) and the Euro- 



5 



GRAPHENE project CONGRAN. 



1 K. S. Novoselov et al, Science 306, 666 (2004). 

2 A. H. Castro Neto et al, Rev. Mod. Phys. 81, 109 (2009). 

3 F. Guinea et al, Nature Physics 6, 30 (2010). 

4 N. Levy et al, Science 329, 544 (2010). 

5 Vitor M. Pereira et al, Phys. Rev. B 80, 045401 (2009). 

6 Z. Ni et al, ACS Nano 3, 483 (2009). 

7 T. M. G. Mohiuddin et al, Phys. Rev. B 79, 205433 (2009). 

8 F. Guinea et al, Nature Physics 6, 30 (2010). 

9 H. Castro Neto et al, Rev. Mod. Phys. 81, 109 (2009). 

10 F. Guinea et al, Phys. Rev. B 77, 075422 (2008). 

11 K-J Kim et al, Phys. Rev. B, 84, 0814401(R) (2011). 

12 S. Philip et al, Appl. Phys. Lett. 94, 032101 (2009). 

13 H. Tomoriet et al, Applied Physics Express 4, 075102 
(2011). 

14 D. W. Brenner et al, J. Phys.: Condens. Matter 14, 783 
(2002). 

15 Zhun-Yong Ong and Eric Pop, Phys. Rev. B 81, 155408 
(2010) 

16 M. Neek-Amal et al. Phys. Rev. E, 82 051605 (2010). 

17 M. Neek-Amal and F. M. Peeters, Phys. Rev. B 81, 235421 



(2010). 

18 A. I. Zhbanov et al, ACS Nano, 4, 5937 (2010). 

19 N. Chandra et al, Phys. Rev. B. 69, 094101 (2004). 

20 Q. X. Pei et al, Carbon 48, 898 (2010). 

21 H. Rafii-Tabar, Physics Reports 390, 235 (2004). 

22 J. W. Klos et al, Phys. Rev. B 80, 245432 (2009). 

23 S. Costamagna, O. Hernandez, and A. Dobry, Phys. Rev. 
B 81, 115421 (2010). 

24 M. J. Schmidt and D. Loss, Phys. Rev. B 81, 165439 

(2010) . 

25 A. Weisse et al, Rev. Mod. Phys. 78, 275 (2006). 

26 G. Schubert and H. Fehske, Phys. Rev. Lett. 108, 066402 
(2012). 

27 L. Covaci, F. M. Peeters, and M. Berciu, Phys. Rev. Lett. 
105, 167006 (2010). 

28 L. Covaci and F.M. Peeters, Phys. Rev. B 84, 24140RR) 

(2011) . 

29 A movie of the LDOS maps for different energies is pro- 
vided in the Supplementary Material. 



